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Abstract 

There are 880 magic squares of size 4 by 4, and 275,305,224 of size 
5 by 5. It seems very difficult if not impossible to count exactly the 
number of higher order magic squares. We propose a method to esti- 
mate these numbers by Monte Carlo simulating magic squares at finite 
temperature. One is led to perform low temperature simulations of a 
system with many ground states that are separated by energy barri- 
ers. The Parallel Tempering Monte Carlo method turns out to be of 
great help here. Our estimate for the number of 6 by 6 magic squares 
is (0.17745 ± 0.00016) x 10 20 . 



1 Introduction 



Magic squares involve using all the numbers 1, 2, 3, . . . , n 2 to fill the squares 
of an n x n board so that each row, each column, and both main diagonals 
sum up to the same number. Interesting information about magic squares is 
collected at the WEB site of M. Suzuki 0. 

We here address the question of the number of magic squares N(ri) of a 
given order n. In the following, this number should always be understood as 
the total number of magic squares divided by 8, thus considering as equiv- 
alent those squares which can be obtained from each other by the obvious 
reflection and rotation symmetries. It has been known since long that there 
are 880 magic squares of order 4. N(5i) was estimated by L. Candy in his 
Construction, Classification and Census of Magic Squares of Order Five, pri- 
vately published in 1938. Candy arrived at a total of 13,288,952. The exact 
number was determined by Richard Schroeppel in 1973, using a computer 
backtracking program, see [0. His result 275,305,224 shows that Candy's 
estimate was low by a wide margin. It seems very difficult to exactly deter- 
mine N(n) for n > 5. However, it is possible to obtain statistical estimates 
with good precision. In this paper we shall describe a Monte Carlo method 
for this purpose. As a demonstration, we apply it to the cases n — 4, 5, 
and 6. Our method is, however, by no means restricted to theses special 
cases, and could well be used for higher n, and also for all kinds of variants 
of magic squares, like pan-magic squares, squares filled with primes only, or 
magic cubes. 

2 Magic Squares at Finite Temperature 

We consider magic squares as the zero temperature configurations of a sta- 
tistical system with partition function 

Z{p) = Y exp[-0E(C)} , (1) 
c 

where the sum is over all possibilities to fill the square with the numbers 
n — 1,2, ... , n 2 . (3 is proportional to the inverse temperature. We define the 
energy of a configuration by 

E(C) = Y (S e - M) 2 + Y, (Sr~ M) 2 + Y (Sa - M) 2 , (2) 

columns c rows r diagonals d 
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where S c , S r , and Sj are the sums of the columns, rows, and main diagonals, 
respectively. M = n(n 2 + l)/2 is called the magic constant. Obviously 
E(C) > 0. Magic squares have zero energy. The task of counting the number 
of magic squares N(n) is thus equivalent to counting the number of states of 
minimal energy, i.e. determining the zero temperature entropy. 

Monte Carlo methods allow to estimate expectation values of functions 
A of the configurations C, 

(A) = WeM-PE(C)}A(C). (3) 
z c 

Now observe that 

N(n) = I lim Z((3) . (4) 
Z(/3) is not an expectation value. However, since Z(0) = (n 2 )!, we have 

N{n) = I Z(0) lim (exp[-/3E(C)]) /3=0 , (5) 

where the subscript indicates that the expectation value has to be taken at 
infinite temperature here. Eq. (pj) is still not practical for calculations, since 
for large (3 the measured quantity exp[— (3E(C)] fluctuates over many orders 
of magnitude, thus leading to very large statistical errors of the the Monte 
Carlo estimate. 

Let us therefore consider a collection of j3- values = (3\ < j3 2 < ■ ■ ■ < Pm- 
Then 

71 ' A = (e-^-**)* , (6) 



Zifli) 
so that 

||| = <e-W"-A>*> A ( C -W.-A)^) A . . . ^-^e^ {e -W-M* )flm . 

(7) 

If the /^-differences in the measured quantities are not too big, this representa- 
tions offers a way to compute Z(j3) for large j3 and thus an approximation for 
N(n). We remark that Z(/3) is strictly monotonously decreasing. The finite 
/9-value therefore yields an upper bound on the number of ground states. 
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3 Parallel Tempering Monte Carlo 



A valid Monte Carlo algorithm to estimate any of the expectation values 
occurring in eq. (0) can be built using the Metropolis procedure: Propose to 
exchange the positions of two entries in the square, determine the correspond- 
ing energy change AE, and implement the modification of configuration with 
probability 

p = rninfl, e~ pAE ) . (8) 

Except for rather small (3, the acceptance rates become prohibitively small 
if one just exchanges randomly selected entries. We therefore restrict the 
update proposal to the transpositions of 1 with 2, then 2 with 3, . . ., n 2 — 1 
with n 2 . Such moves are also useful in Simulated Annealing procedures 
designed to search for magic squares with very large n, by minimizing E(C) 
or a similar cost function. 

For larger /3-values, naive Monte Carlo simulations run into a problem: 
The different areas of low energy are separated by high barriers. In order to 
sample all the low energy contributions with the right weight one has to pen- 
etrate and tunnel through these barries. Consequently, very long simulation 
times are needed. 

The situation can be much improved by using the Parallel Tempering 
or Exchange Monte Carlo method, see, e.g. and further references cited 
in [0J. It amounts to simulate the joint ensemble of all the (independent) 
systems with inverse temperatures in parallel. The partition function of 
this system is 

^joint = E E • • • E exp[-/3 1 E(C 1 ) - (3 2 E{C 2 ) - ... - (3 m E{C m )} . (9) 

C\ C*2 Cm 

In addition to updating independently the configurations Cj, one includes 
exchanges of configurations, usually of adjacent f3- values. A proposal of such 
a change is again accepted using a Metropolis procedure. E.g., configurations 
Cj and Cj + i are exchanged with probability 

Pm+1 = rninfl, e -^ + ^)(z(C)-E(c. +1 y }] (1Q) 

The exchange of configurations over the temperature range strongly speeds 
up the Monte Carlo process at the lower temperatures. Numerical experience 
shows that, in order for the procedure to be efficient, the acceptance rates for 
the configuration exchanges should not be very much smaller than one half. 
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Furthermore, having too many systems might also hamper rapid exchange 
of information from higher to lower temperatures and vice versa. 

4 Monte Carlo Results 

We simulated squares with n = 4, 5 and 6, using always a set of m = 20 
/3-values. The results are summarized in the tables [3], @, and |3|. The largest 
/?-value was chosen such that the acceptance rate ui m = uj{(3 m ) was around 
one percent. The intermediate /5-values were chosen such that the exchange 
rates in the tempering cycle are roughly of order one half or bigger. Each 
tempering cycle consisted in performing one Metropolis updating sweep for 
each of the 20 configurations and then attempting to exchange each of the 
adjacent /3j-pairs. We made 3.25- 10 7 for n = 4, and 10 s such cycles for n = 5 
and n = 6, respectively. This required approximately a total of 12 days on 
a 166 MHz Pentium PC. The code was not optimized yet with respect to 
run-time behaviour. 

The Pi, the acceptance rates u>i at Pi, and the exchange rates oJ^ i+ i are 
given in columns 2, 3, and 4 of the tables. The energy expectation value 
estimates are given in column 5. The last columns of the tables give the 
ratios of partition functions Z (P i+ i) / Z (Pi) , where Zi = Z(Pi). The bottom 
parts present Z(P)/Z(p 20 ) for three extra /3- values much larger than P m . 
Having three /3's of increasing size allows us to check for convergence of 
the N(n). The errors of these estimates were obtained by generating 50 
synthetic data sets of the Zi + \/Zi, scattering them around the measured 
values according to a Gaussian distribution with variances given by the error 
bars of the simulation results. 

Both for n = 4 and n = 5, our estimates for N(n) stabilize reasonably 
with increasing P and agree with the exactly known results. Our estimate 
for JV(6) is 0.17745(16) ■ 10 20 . 

5 Conclusions 

The method proposed provides reliable estimates for the numbers of magic 
squares. Of course, there remain many possibilities to improve on the simu- 
lations (besides going to higher statistics). E.g., one could play around with 
m, the choices of the Pi and also with the frequency of configuration exchange 
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attemps. 

It would be interesting to complement the present approach with, e.g., 
analytical methods. High temperature expansions seem feasable. For exam- 
ple, the energy at f3 = can be fairly easily evaluated exactly, and is given by 
(E)q = n 2 (n 4 — l)/6. One can convince oneself that the higher moments of 
the energy at infinite temperature can all be expressed in closed form, most 
likely as polynomials in n. 

A very short run (10 6 cycles) for n = 7 yields N(7) = 0.3760(52) • 10 35 . It 
is an interesting question whether there is some simple behaviour of N(n). 
Finally, it could be worthwile to study much larger n to look out whether 
the magic squares at finite temperature have also interesting thermodynamic 
properties like phase transitions. 
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n 


= 4 


3.25 • 


LO 7 tempering cycles 


i 


Pi 








Zi+i/Zi 


1 


0.0000 


1.000 


0.868 


680.048(56) 


0.520412(29) 


2 


0.0010 


0.990 


0.922 


626.364(50) 


0.693277(20) 


3 


0.0016 


0.984 


0.861 


594.757(51) 


0.536070(29) 


4 


0.0027 


0.973 


0.800 


539.585(49) 


0.427069(34) 


5 


0.0044 


0.957 


0.713 


463.419(45) 


0.315796(35) 


6 


0.0072 


0.934 


0.631 


365.548(32) 


0.234237(28) 


7 


0.0119 


0.902 


0.582 


262.108(28) 


0.196354(33) 


8 


0.0195 


0.862 


0.552 


176.937(21) 


0.172897(33) 


9 


0.0319 


0.814 


0.529 


114.919(11) 


0.157632(26) 


10 


0.0523 


0.754 


0.516 


72.5736(73) 


0.150015(30) 


n 


0.0858 


0.681 


r\ r AO 

0.508 


44.9679(50) 


0.148085(27) 


1 o 


n 1 /iA7 
U.14U ( 


U.594 


U.5UU 


27.3797(28) 


0.151383(29) 


1 Q 


U.zoUo 


U.4y 1 


n /io/i 


16.2867(17) 


0.163141(32) 


14 


0.3786 


0.374 


0.495 


9.3625(11) 


0.188879(35) 


15 


0.6208 


0.248 


0.485 


5.08151(74) 


0.264513(67) 


16 


1.0000 


0.118 


0.704 


2.29517(72) 


0.60598(12) 


17 


1.3000 


0.057 


0.813 


1.15538(60) 


0.78174(11) 


18 


1.6000 


0.026 


0.925 


0.55536(40) 


0.91600(60) 


19 


1.8000 


0.015 


0.951 


0.33992(26) 


0.94757(40) 


20 


2.0000 


0.009 




0.20941(20) 






P 






z(P)/z 20 


7V(/3)/10 3 




5.0 






0.912727(81) 


0.87968(57) 




8.0 






0.912572(81) 


0.87953(58) 




11.0 






0.912571(81) 


0.87953(58) 




exact 








0.880 



Table 1: Monte Carlo results for n = 4. 
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n = 5 


10 8 tempering cycles 


% 


A 








Zi+i/Zi 


1 


0.0000 


1.000 


0.561 


2600.09(17) 


0.10558(15) 


2 


0.0010 


0.982 


0.730 


1585.362(97) 


0.29389(21) 


3 


0.0017 


0.971 


0.633 


1188.134(67) 


0.19578(19) 


4 


0.0029 


0.955 


0.551 


816.665(45) 


0.13749(16) 


5 


0.0049 


0.934 


0.495 


527.305(33) 


0.10641(13) 


6 


0.0084 


0.907 


0.462 


327.343(19) 


0.09006(13) 


7 


0.0142 


0.873 


0.442 


198.631(11) 


0.08139(10) 


8 


0.0241 


0.831 


0.431 


118.9178(64) 


0.07663(10) 


9 


r\ r\ A -1 r\ 

0.0410 


0.777 


0.424 


70.6331(35) 


0.07403(10) 


10 


0.0698 


0.710 


0.421 


41.7659(21) 


0.0726609(97) 


11 


0.1186 


0.626 


0.419 


24.6289(12) 


0.0719532(95) 




0./016 


O.OZO 


0.418 


14.48913(70) 


0.0716943(94) 


1 Q 


n Q/107 

U.o4z / 


U.4Uy 


0.418 


8.50502(44) 


0.0718681(90) 


14 


0.5826 


0.285 


0.415 


4.90719(32) 


0.073181(11) 


15 


0.9905 


0.165 


0.347 


4.90719(32) 


0.093419(26) 


16 


1.6838 


0.058 


0.810 


2.21324(49) 


0.659590(71) 


17 


1.9000 


0.039 


0.779 


1.65545(52) 


0.66774(10) 


18 


2.2000 


0.021 


0.881 


1.07235(50) 


0.830839(79) 


19 


2.4000 


0.014 


0.903 


0.79392(47) 


0.871964(76) 


20 


2.6000 


0.009 




0.58657(42) 






P 






z(P)/z 20 


N(f3) /10 9 




5.6000 






0.66293(25) 


0.27914(19) 




8.6000 






0.65469(26) 


0.27577(19) 




11.6000 






0.65429(26) 


0.27550(19) 




exact 








0.275305204 



Table 2: Monte Carlo results for n = 5. 
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n = 6 


10 s 


tempering cycles 


i 


0i 






(E)r 




1 


0.0000 


1.000 


0.513 


7768 23(61) 


070860(14) 

V 7 • V J 1 V 7 v7 V7 V 7 \ -L J- j 


2 


0.0004 


0.988 


0.638 


5615 65(36) 


168908(18) 


3 


0.0008 


0.979 


0.502 


4352 01(27) 


086142(13) 


4 


0.0014 


0.966 


0.403 


2976 04(19) 


0494970(93) 

V 7 • V 7 J- '.7 -L '.7 1 V 7 \ ' 7 v J 


5 


0.0027 


0.948 


0.345 


1827.89(12) 


0.0342720(72) 


6 


0.0052 


0.924 


0.314 


1046.796(62) 


0.027587(58) 


7 


0.0099 


0.892 


0.298 


576.290(33) 


0.0244279(51) 


8 


0.0188 


0.849 


0.289 


310.555(17) 


0.0228945(51) 


9 


0.0358 


0.791 


0.285 


165.4935(79) 


0.0221136(43) 


10 


0.0679 


0.713 


0.283 


87.6788(38) 


0.0216969(42) 


11 


0.1291 


0.612 


0.282 


46.3095(20) 


0.0214833(42) 


12 


0.2452 


0.486 


0.281 


24.4166(11) 


0.0213753(44) 


13 


0.4660 


0.340 


0.282 


12.86353(59) 


0.0213645(47) 


14 


0.8853 


0.193 


0.235 


6.74920(32) 


0.0262622(88) 


15 


1.6821 


0.065 


0.790 


2.91925(39) 


0.570565(54) 


16 


1.9000 


0.045 


0.824 


2.25025(44) 


0.672457(64) 


17 


2.1000 


0.031 


0.844 


1.73410(48) 


0.738535(73) 


18 


2.3000 


0.021 


0.930 


1.31250(49) 


0.884956(43) 


19 


2.4000 


0.018 


0.877 


1.13603(48) 


0.821399(76) 


20 


2.6000 


0.012 




0.84515(46) 






P 






z(p)/z 20 


N((3)/10 20 




6.6000 






0.55845(25) 


0.17842(16) 




10.6000 






0.55542(25) 


0.17745(16) 




14.6000 






0.55536(25) 


0.17744(16) 



Table 3: Monte Carlo results for n = 6. 
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